Theory of optical axion electrodynamics and application to the Kerr effect in topological antiferromagnets

Emergent axion electrodynamics in magneto-electric media is expected to provide novel ways to detect and control material properties with electromagnetic fields. However, despite being studied intensively for over a decade, its theoretical understanding remains mostly confined to the static limit. Here, we introduce a theory of axion electrodynamics at general frequencies. We define a proper optical axion magneto-electric coupling through its relation to optical surface Hall conductivity and provide ways to calculate it in lattice systems. By employing our formulas, we show that axion electrodynamics can lead to a significant Kerr effect in thin-film antiferromagnets at wavelengths that are seemingly too long to resolve the spatial modulation of magnetism. We identify the wavelength scale above which the Kerr effect is suppressed. Our theory is particularly relevant to materials like MnBi2Te4, a topological antiferromagnet whose magneto-electric response is shown here to be dominated by the axion contribution even at optical frequencies.


Junyeong Ahn 1 , Su-Yang Xu 2 & Ashvin Vishwanath 1
Emergent axion electrodynamics in magneto-electric media is expected to provide novel ways to detect and control material properties with electromagnetic fields. However, despite being studied intensively for over a decade, its theoretical understanding remains mostly confined to the static limit. Here, we introduce a theory of axion electrodynamics at general frequencies. We define a proper optical axion magneto-electric coupling through its relation to optical surface Hall conductivity and provide ways to calculate it in lattice systems. By employing our formulas, we show that axion electrodynamics can lead to a significant Kerr effect in thin-film antiferromagnets at wavelengths that are seemingly too long to resolve the spatial modulation of magnetism. We identify the wavelength scale above which the Kerr effect is suppressed. Our theory is particularly relevant to materials like MnBi 2 Te 4 , a topological antiferromagnet whose magneto-electric response is shown here to be dominated by the axion contribution even at optical frequencies.
In a medium where spatial inversion and time reversal symmetries are both broken, magnetic and electric fields couple in a way that is fundamentally different from the electromagnetic induction described by Maxwell's equations in vacuum, thus leading to exotic electrodynamics. A topic of particular interest in magneto-electric coupling phenomena recently is axion electrodynamics. Axion electrodynamics is theoretically proposed to be realized in a class of topological materials called axion insulators 1-6 for sufficiently slowly oscillating electromagnetic fields. In particular, the static axion magneto-electric coupling is quantized by a multiple of the fundamental constant e 2 / 2h [7][8][9] , which originates from the half-quantized Hall conductance of the topological surface states. While such a quantized magneto-electric coupling has not yet been directly observed, experimental progress has been made to observe its consequences in the static limit [10][11][12][13][14][15][16] .
In contrast to the static limit, the emergent axion electrodynamics at generic optical frequencies remains largely unexplored theoretically. Formulating a theory of optical axion electrodynamics has to overcome the challenge that the axion coupling is ill-defined in periodic three-dimensional systems, making it hard to calculate and understand. The static axion angle can be calculated by the Chern-Simons integral of the non-abelian Berry connection 8,17,18 at the cost of being gauge dependent. Despite its gauge dependence, the Chern-Simons integral is well-defined as an angular variable taking a value between 0 and 2π, because a gauge transformation changes its value only by a multiple of 2π. This is a manifestation that the static surface Hall conductance changes by a multiple of e 2 /h under deformations. However, at optical frequencies, a similar approach does not seem feasible. This is because the optical surface Hall conductivity is frequency dependent, and thus a generic surface deformation changes the surface Hall response by a non-quantized amount. Owing to this difficulty, theoretical understanding of the optical axion electrodynamics has remained elusive to date 19 .
In this paper, we make two important steps toward the complete formulation of the optical axion electrodynamics. First, we show that a proper definition of the optical magneto-electric coupling allows us to calculate the optical axion angle in a fully gauge-independent way in a system with finite thickness. Although the optical axion electrodynamics is a part of the linear-response optical magnetoelectric effects, it is distinguished from the other contributions. Nonaxionic magneto-electric effects are described within the well-established theory of gyrotropic birefringence and natural optical activity, which are based on the bulk current response [19][20][21] . On the other hand, the optical axion electrodynamics is a surface phenomenon and thus is not captured by those theories 19 . We therefore define the well-define optical axion angle by analyzing the surface response. Using the well-defined axion angle, we can understand the optical axion electrodynamics of thin films and also that of bulk crystals by increasing the thickness. Second, in systems with periodic boundary conditions along all three directions, we find the optical axion angle at high optical frequencies can be estimated by the optical layer Hall conductivity that we define in the main text. These findings present advantages in numerical calculations as well as conceptual advances.
Our development of a theory of optical axion electrodynamics opens the door to understanding novel axion-induced optical phenomena. Here we provide a concrete example. Magneto-optic effects make electromagnetic waves powerful probes of the magnetic structure in materials. Conventionally, both Kerr and Faraday effects are attributed to the optical Hall effect and thus commonly perceived as probes of net magnetization 22,23 . Although not well known, previous theoretical studies proposed that the magneto-electric effect can also lead to the Kerr effect, providing a novel way to probe fully compensated antiferromagnets whose symmetries strictly prohibit the Hall effect and net magnetization 20,[24][25][26][27][28] . Still, however, there are two aspects that need further investigation. First, axion electrodynamic contribution to the Kerr effect has not been well understood. Second, since previous studies have focused on three-dimensional bulk systems, magneto-optic Kerr effects in quasi two-dimensional antiferromagnets remain elusive. We present theoretical analysis revealing the precise conditions for realizing the optical axion electrodynamic Kerr effect as well as quantitative numerical analysis allowed by our gaugeinvariant formulas.
Our results apply broadly to both topological and non-topological media because optical magneto-electric effects, as non-quantized phenomena, are not sensitive to the topological nature of the ground state. Therefore, our work provides a theoretical basis for the detection and manipulation of antiferromagnetism in a large class of materials, thus having potential for wide applications to antiferromagnetic spintronics and the study of magnetic structure in quantum materials Meanwhile, topological antiferromagnets need special attention as they are ideal platforms for optical axion electrodynamics. To understand this, we note two aspects. First, it is desired to have antiferromagnets having bulk symmetry that reverses an odd number of spacetime coordinates and is broken on the surfaces (e.g., spatial inversion symmetry), because such a symmetry suppresses bulk magneto-electric effects but not the axion magneto-electric effect. This condition is similar to the requirement of the quantization of the static axion angle but additionally requires that the quantizingsymmetry is broken at the surface to allow for a nonzero value. Second, spatially spreading of electronic quantum states is needed, in order to show a response distinguished from that of decoupled layered or Mott antiferromagnets. This condition is again satisfied in topological antiferromagnets. We thus apply our theory to a model of MnBi 2 Te 4 , which is the only experimentally realized axion topological antiferromagnet to date, and show that the Kerr effect in this material is significant.

Optical magneto-electric coupling
Motivated by the equivalence between the surface Hall conductivity and the axion angle in the static limit, we study the surface current response to define the optical axion angle. We consider currents generated by electromagnetic multipole moments. Electromagnetic responses from multipole moments are smaller for higher order moments 21 : electric dipole P i ≫ electric quadrupole Q ij and magnetic dipole M i ≫ higher orders. Here, we consider only up to the electric quadrupole-magnetic dipole order, giving the leading-order magnetoelectric effect. The bulk current density is related to the multipole moments by: While electric quadrupole and magnetic dipole moments do not generate macroscopic currents in macroscopically homogeneous lattice systems, they generate currents on the system boundary where the material property changes spatially. The induced multipole moments have the following form in the frequency domain 21 : where χ ij = χ ji and χ 0 ij = À χ 0 ji are the electric susceptibility tensors, G ij and G 0 ij are magneto-electric coupling, and a 0 ijk = a 0 ikj and a ijk = a ikj are electric quadrupolar susceptibility tensors. Here, we are interested in the magnetic magneto-electric effects described by G ij and a 0 ijk , which occur only when time reversal symmetry is broken while being compatible with spacetime inversion PT symmetry. Therefore, we assume PT symmetry to neglect the complications arising from the bulk Hall conductivity and natural optical activity, both of which are excluded in this symmetry setting, i.e., χ 0 ij = 0 and G 0 ij = a ijk = 0 (Table 1). This assumption does not affect our key results (see Methods for discussions without PT symmetry). As we show below, magnetoelectric effects occur in combination with electric quadrupole responses at optical frequencies, requiring the consideration of the combination of G ij and a 0 ijk . Let us consider the surface at z = 0 of a three-dimensional homogeneous material with the outward normal directionẑ (Fig. 1). The surface current density j s = R d s =2 Àd s =2 dzJ, where d s is the characteristic thickness of the interface where response functions change rapidly as functions of z, is related to the multipole moments through where δf = f( − d s /2) − f(d s /2) and f = ½f ðÀd s =2Þ + f ðd s =2Þ=2 are the difference and average of the material property f(z) across the interface, and σ ij and σ ijk = i(ϵ jkl T il + ϵ ikl T jl ) + ωS ijk are the bulk conductivity tensors defined by We neglect the δE term because it is a bulk response whose contribution to the interface vanishes as d s → 0. The form of the surface current density in Eq. (2) suggests defining surface-sensitive magnetoelectric coupling by where the superscript (z) explicitly shows that these quantities depend on the choice of the surface normal direction [(z) and (−z) are the same, though]. The zz component of this surface-sensitive magneto-electric coupling is defined from the response of the surface charge. Using ρ = À ∇ Á P + 1 2 ∂ i ∂ j Q ij + . . ., we obtain ρ s = G zi B i + ω 2 ða 0 yzx À a 0 xzy ÞB z + ω 2 a 0 zzy B x À ω 2 a 0 zzx B y + . . . where the ellipsis includes the electric dipole term χ zi E i , the symmetric (∂ j E k + ∂ k E j ) terms, and higher-order multipole terms. From this, we define as well as G ðzÞ zx = G zx + ωa 0 zzy =2 and G ðzÞ zy = G zy À ωa 0 zzx =2. While various components of the magneto-electric coupling are related to the surface optical conductivity, there is a unique component that manifests itself only at the surface: the axion angle. We define the optical axion angle by the trace part of G ðzÞ ij .
θ ðzÞ ðωÞ π 2h e 2 1 3 To see that this only manifests itself at the surface rather than the bulk, we note that the bulk current response by the magneto-electric and electric quadrupole susceptibilities is determined only through T ij and S ijk . Because of this, previous studies focusing on bulk magnetoelectric effect did not capture optical axion electrodynamics 19 Since T ij is traceless, it does not contribute to the trace of G ðzÞ ij . Note that G ðzÞ ij ðωÞ transforms as a tensor under magnetic layer group actions but not under the full three-dimensional magnetic space group actions.

Magneto-electric coupling with open boundaries in one direction
Defining the combination G ðzÞ ij of G ij and a 0 ijk has an advantage in practical calculations as well as in the formulation. As G ðzÞ ij characterizes the surface current response which is measurable, it admits a gaugeinvariant form in three-dimensional lattice systems with open boundaries along one direction, or in quasi-two-dimensional systems. This nice property is absent in the diagonal magneto-electric coupling G ii , making it hard to calculate.
The bare magneto-electric coupling G ij = ∂ Pi /∂B j and a 0 ijk = ∂Q jk =∂E i at zero temperature have the following form according to linear response theory (see Methods).
where V is the volume of the system, n, m are indices for energy eigenstates with eigenvalue ℏω n , f nm = f n − f m is the difference between the Fermi-Dirac distribution of the n and m states.P i = À er i =V , magnetic dipole, and electric quadrupole density operators, respectively, wherer andv are the position and velocity operators of electrons, andm s is the spin magnetic moment operator. Equation (7) can be calculated for molecular systems 21,29 , meaning finite systems with open boundary conditions along all directions, by using the realspace representation ofr and the relationv = À i_ À1 ½r,Ĥ, whereĤ is the Hamiltonian of the system. In periodic systems, however, G ij and a 0 ijk are not well defined separately because the position operator is not well-defined because the position is not uniquely defined. This manifests through the momentum-space representation of the position operator hψ mk 0 |r|ψ nk i = À δ mn i∂ k δ k 0 k + δ k 0 k hu mk |i∂ k |u nk i, whose diagonal matrix element is not well defined because of ∂ k δ k 0 k . On the other hand, the diagonal matrix elements of the position operator do not appear in the response functions that have a well defined physical meaning in periodic systems. T ij and S ijk are such examples that characterizes the bulk current response 19 . Since G (z) characterizes the surface response, one can expect that it is well defined in quasi-two-dimensional periodic systems.
After doing some algebra that we relegate to Methods, we can write G ðzÞ ii s in two-dimensional momentum space as where r i mn = hψ m |r i |ψ n i and v i mn = hψ m |v i |ψ n i are matrix elements of the position and velocity operators, |ψ nðk x ,k y Þ i is the Bloch state, the subscripts n,m and p are band indices. While the position operator matrix element is not well defined in momentum space 30 , the diagonal matrix elements ofr x andr y do not appear in Eq. (8), so that the surfacesensitive magneto-electric coupling is well defined in two-dimensional momentum space. Combining equations in Eqs. (5) and (8), we can obtain the optical axion angle. In simple tight-binding models wherer z commutes withv y , G ðzÞ xx ð0Þ reduces to the expression of G xx (0) derived by Liu and Wang 31 . To our knowledge, our formulas represent the first expressions for calculating the diagonal components of the optical magneto-electric coupling in crystalline (periodic) systems.
The form of G ðzÞ xx and G ðzÞ yy suggests that their orbital magnetoelectric part may be interpreted as a real-space dipole of the Berry curvature. This idea works exactly for two-band systems in the limit of decoupled layer systems, where the matrix element part can be written as the product of the Berry curvature F xy = À 2Im½r x 12 r y 21 and the z component of the position eigenvalues r z 11 or r z 22 . We apply this idea below to the case where the z direction is also periodic.

Intra-cell magneto-electric coupling and optical layer Hall conductivity
In fully periodic lattice systems where the z direction is also periodic, it is hard to calculate the full orbital part of G ðzÞ ii ðωÞ becauser z is not well defined in momentum space. Nevertheless, we can still define and calculate the magneto-electric coupling of the unit cell, which we call the intra-cell magneto-electric coupling. Note, this treatment will be necessarily approximate, in contrast to our previous discussion of the slab geometry, but provides an physical understanding for the results. For example, the use of the intra-cell magneto-electric coupling makes the relation concrete between the axion angle and the antiferromagnetism in inversion symmetric systems. Ultimately we must compare the results between this approach and the previous slab calculation as we do in another section below.
Let us begin by giving a physical intuition that a nonzero axion angle is natural in a fully compensated antiferromagnet 32 . To see this, recall the analogy between electric polarization in 1D and the axion theta angle θ in 3D. The former can be defined in a system free of net charge, by stacking alternate positive and negative charges. Similarly, if we stack alternate planes with Hall conductance ± g H xy such that the net Hall conductance vanishes, the axion magneto-electric coupling becomes well defined and is the analog of electric polarization of the alternating planes. In fact, with an applied magnetic field perpendicular to the planes, the induced polarization from having alternating charges in the ± g H xy layers, leads to a finite electric polarization which is readily calculated as g H xy δa z =a z , where δa z is spacing between alternating antiferromagnetic planes versus the vertical size of the unit cell a z . The charge in each flux quantum area is g H (h/e). For an antiferromagnet where spacing between planes = 1/2a z , this is just θ = 2πg H xy =2. To present a more detailed analysis, let us decompose the position operator into the intra-cell polarization and unit cell position parts by r z mn = A z mn + R z mn . Namely, where |w αR i is the Wannier state with the collective index α for spin and orbital (cf. n and m are band indices), |ψ αk i = N À1=2 P R e ikÁR |w αR i, and A z βα ðkÞ = P R hw β0 |r z |w αR ie ikÁR . The second term in Eq. (9) vanishes for n ≠ m as one can see by writing it as Àiδ mn ∂ k z δ k,k 0 . This decomposition is independent of choosing the basis α within the unit cell, where each unit cell is labeled by a given lattice vector R. However, the decomposition depends on the choice of the unit cell, which we discuss below. The intra-cell polarization term, the first term in Eq. (9), defines the magneto-electric coupling of the unit cell. This intra-cell magnetoelectric coupling depends on the choice of the unit cell. In our case, choosing a unit cell corresponds to fixing the value of the Wannier centers hw α0 |r z |w α0 i, which is ambiguous by respective lattice translations of the Wannier states |w α0 . A physical interpretation of this multi-valuedness is similar to that of electric polarization 32 : The magneto-electric coupling depends on how we open the boundary, and there exists a preferred choice of the unit cell for each boundary condition ( Fig. 2(a)).
The unit-cell position term has the form P R z R z σ H xy ðR z Þ, where σ H xy ðσ xy À σ yx Þ=2 is the Hall conductivity, and σ H xy ðR z Þ = Àe 2 _ À1 V À1 P n≠m,k x ,k y f nm ω mn =ðω mn À ωÞ À1 Im½r x nm r y mn P R z nn , andP This term is also multivalued in periodic systems because the unit cell position R is not uniquely defined in periodic systems. This part, however, does not contribute to the magneto-electric coupling when the Hall conductivity of the unit layer vanishes, which is the case in PT-symmetric systems. When the boundary is introduced, the Hall conductivity of the surface unit layer can be nonzero even when the bulk Hall conductivity vanishes ( Fig. 2(b)). This change of the R z term is the main source of the difference in the magneto-electric couplings between finite-size systems and periodic lattice systems. This effect is significant especially in the static limit of axion insulators, where the emergent Dirac surface states modify the lowfrequency surface Hall conductivity in the order of e 2 /2h.
As we derive below, in inversion-symmetric even-layer antiferromagnets, the intra-cell contribution to the magneto-electric coupling is simplified to the optical layer Hall conductivity. Namely, the optical axion angle is where a z is the vertical lattice constant, and σ xy is a three-dimensional conductivity taking the unit of conductance (unit of e 2 /h) divided by the length. As for the definition of layers, we note that there are two inversion-invariant values of the vertical displacement z within a z . We refer to the two quasi-two-dimensional bipartite regions centered at those two invariant z values as layers ( Fig. 2(a)). The layer Hall conductivity is then defined as σ LH xy = ðσ H,u xy À σ H,d xy Þ=2 from the bulk optical Hall conductivity projected to the single layer l = u, d: where P l n 0 n ðkÞ = P α,R;r z α,R 2l layers hψ nk |w αR ihw αR |ψ n 0 k i, |w αR is are inversion-symmetric Wannier states, and r z α,R = hw αR |r z |w αR i is the Wannier center.
To obtain Eq. (10), note that the unit cell for even-layer systems is symmetric under the combination of inversion and a lattice translation of either the u layer by −a z (or the l layer by + a z ) ( Fig. 2(a)). If we focus on the intra-cell part (the A part), the combined symmetry gives a constraint θ ðzÞ A ðωÞ = À θ ðzÞ A ðωÞ À a z σ H,d xy ðωÞ, where the last term is due to the transformation property of the axion angle δθ ðzÞ ðωÞ = d z σ H xy ðωÞ under the translation by a vector d. As we consider antiferromagnets with zero net Hall conductivity, such that σ H,d xy = À σ H,u xy = À σ LH xy , we arrive at Eq. (10). Another useful way of understanding this result is to think of the electric polarization generated by applying a uniform magnetic field, and then the displacement of layers of one sign of the Hall effect obviously contributes to change in electric polarization, as we explained at the beginning of this section.

Magneto-optic effects in fully compensated antiferromagnets
The optical axion angle manifests directly through the magneto-optic Kerr effect. To gain some intuition for this, recall that the change in axion angle at the surface leads to a surface Hall effect σ s xy = ðδθ=2πÞe 2 =h. Clearly, such a surface Hall conductance will lead to a Kerr effect. We present a more systematic analysis below.
Let us consider the reflection at the single interface between two media 1 and 2 (Fig. 3). For simplicity, we assume that both media have M x T, C 3z , and PT symmetries, which are shared by magnetoelectric materials Cr 2 O 3 and MnBi 2 Te 4 . We also assume normal incidence (i.e., light is incident along Àẑ). Then there is no birefringence because of the symmetry we require, so the propagation of light within each medium is then characterized by a single complex-valued refractive index n μ , where μ = 1, 2 labels the two media (See Eq. (6)).
The electric field in medium 1 consists of incident and reflected fields E i and E r while that in medium 2 is the transmitted field E t : where r and t = 1 + r are 2 × 2 Jones matrices for reflection and transmission, respectively. E 1 = E 2 at the interface because electric fields are continuous by Faraday's law, because we consider normal incidence such that electric fields are parallel to the interface. The contribution from the surface conductivity is encoded in the boundary condition of the magnetic field at the interface.
Here, we consider only the surface currents induced from the bulk, i.e., . By solving the boundary condition equations, we obtain r xx = ðn 1 À n 2 Þðn 1 + n 2 Þ À ðμ 0 cσ s xy Þ 2 ðμ 0 cσ s xy Þ 2 + ðn 1 + n 2 Þ 2 , r xy = À 2n 1 ðμ 0 cσ s xy Þ ðμ 0 cσ s xy Þ 2 + ðn 1 + n 2 Þ 2 : Magneto-electric coupling appears in the reflective Jones matrix through the surface conductivity, which is a manifestation that the Kerr effects here are surface phenomena. The complex Kerr angle is defined by ϕ K = φ K + iη K = tan À1 ðr xy =r xx Þ. Its real part measures the rotation of the light  polarization plane while its imaginary part measures the circular dichroism, i.e., the intensity imbalance between the reflected left and right circularly polarized light. The Kerr angle is Oðμ 0 cσ s xy Þ in general, but it can be much enhanced when n 1 ≈ n 2 because of the suppression of r xx .
When the sample thickness is much larger than the wavelength of light, it is enough to suppose that reflection occurs at a single interface. However, when thickness d is comparable to or less than wavelength λ, which is the case particularly relevant for thin films, we need to consider the response from the whole sample including top and bottom surfaces (cf. When the photon energy is smaller than the band gap, double interfaces can be relevant even for d ≫ λ for an insulating medium 14 because then light incident on the top reaches the bottom without attenuation.).
To study this case, we consider three media with refractive index n μ , where μ = 1, 2, 3, as shown in Fig. 3b. We again assume M x T, C 3z , and PT symmetries for each media. The Jones reflection matrix of the sample for light incident from medium 1 is then the infinite sum of multiple reflections.
where we use t T,B = 1 + r T,B and t 0 T = 1 + r 0 T , and ϕ = n 2 ωd/c is the complex-valued phase obtained by the one-way propagation through the sample thickness d. Here, the subscripts T and B indicate the top and bottom of the sample, and the prime indicates the process where the light is incident to the top surface from below (we follow the notation in ref. 13). See Methods for the expression of Jones matrices r T , r 0 T , and r B . Similarly, for transmission, Note that t ≠ 1 + r when ϕ ≠ 0.
Considering the case where only medium 2 is magneto-electric, we set G ðzÞ,1 xx = G ðzÞ,3 xx = 0 and G ðzÞ,2 xx G ðzÞ xx ≠ 0. We first consider |δn| ≪ |ϕ|, where δn ≡ n 3 − n 1 . The Kerr and Faraday rotation angles are then given by tan ϕ K = 2n 1 μ 0 cG ðzÞ xx n 2 2 À n 2 1 + ðμ 0 cG ðzÞ xx Þ 2 + O δn ð Þ, tan ϕ F = À μ 0 cG ðzÞ xx sin ϕ n 2 1 + n 2 2 + ðμ 0 cG ðzÞ xx Þ 2 sin ϕ + 2in 1 n 2 cos ϕ where tan ϕ F = t xy =t xx . A nonzero Faraday rotation requires δn ≠ 0 because PT symmetry needs to be broken 33 . On the other hand, the Kerr rotation can be nonzero with δn = 0 because it is compatible with PT symmetry ( Table 1). The Kerr rotation is independent of ϕ in the leading order of δn when δn is sufficiently small as if the response comes from the top surface only, while it actually comes from multiple reflections between the top and bottom surfaces. In this limit, the way the Kerr effect goes away is highly nontrivial as the thickness goes to zero (i.e., |ϕ| → 0 with |ϕ| ≫ |δn|). The Kerr angle stays constant, but the amplitudes of the reflected signals ultimately vanish. Therefore, a finite Kerr effect can be observed from a thin film when optical equipment is highly sensitive. However, a nontrivial ϕ dependence appears when ϕ is the smallest parameter (|δn| ≫ |ϕ|), where we have which show that both Kerr and Faraday rotation vanish for ϕ = 0. Therefore, a nonzero ϕ is necessary for the Kerr effect as well as the Faraday effect in fully compensated antiferromagnets. Note that, while the Kerr angle in Eq. (17) is independent of ϕ, it applies only when |ϕ| ≫ |δn|.
The crossover between two regimes respectively described by Eqs. (17) and (18) occurs when |ϕ|~|δn|. The corresponding wavelength (n 2 ,T 2 ,S 2 ) (n 1 ,T 1 ,S 1 ) e iφ e iφ e iφ e iφ (n 2 ,T 2 ,S 2 ) (n 1 ,T 1 ,S 1 ) (n 3 ,T 3 ,S 3 ) is much larger than the sample thickness d when |δn| ≪ 1. This is remarkable because one might naively expect that, because the response from the spatially separated top and bottom layers are not resolved when λ ≫ d, the Kerr effect is vanishingly small in that regime and scales linearly with d/λ. However, our analysis shows that a much stricter λ ≫ λ * is required for the suppression of the Kerr angle. To explain how this works, we note again that the Kerr angle and the amplitude of the Kerr rotated signal can behave differently. While the amplitude of the Kerr rotated signal ( ∝ r xy ) is indeed suppressed for λ ≫ d, the amplitude of the non-rotated signal (∝r xx ) is also suppressed in the same limit. Their suppression at large wavelengths compensate each other to keep the ratio ϕ K = r xy /r xx as long as λ is smaller than a larger length scale λ * , above which the Kerr angle ultimately gets suppressed.
In thin film axion insulators, Eq. (18) is typically more relevant at photon energies below the surface gap. The surface gap of experimentally realized axion insulators are about 50 meV 16,34,35 . For example, if we take n 2 = 5 and d = 1 nm, we obtain a very small value |ϕ| ≤ 1.27 × 10 −3 at ℏω ≤ 50 meV, which is typically smaller than |δn|. On the other hand, in the infrared and visible regime where the photon energy is in the order of eV, equation (17) can become relevant. We demonstrate this in the following section.

Model calculations
Let us apply our theory to study the optical axion electrodynamics in MnBi 2 Te 4 . MnBi 2 Te 4 is the only stoichiometric compound that experimentally realizes the intrinsic antiferromagnetic axion insulator 34,35 , which has now become an attractive platform for studying axion magneto-electric effects [10][11][12][36][37][38][39] . As it is a layered antiferromagnet, its few-layer behavior and layer number dependence is also of interest 31,40 .
Here we calculate its magneto-optic properties based on the lowenergy model in refs. 36,41. The goal of our calculations here is to cement the validity of our new theory by providing a concrete model example as well as to understand qualitative features (e.g., the dominance of the axion contribution and the significant Kerr and negligible Faraday effects) of the magneto-optic response in MnBi 2 Te 4 . Our model is expected to quantitatively capture the low-energy properties of the material. On the other hand, at photon energies much larger than the band gap, a precise quantitative calculation of the magnetooptical spectrum will require a model, like the full ab-initio model, that captures all the significant optical transitions involving those states neglected in our model. With this in mind, we consider a nearestneighbor tight-binding Hamiltonian on a layer-stacked triangle latticê where i, j are the site indices, i,j means that the summation is over nearest neighbors, and α, β = 1, …, 4 run over two spin and two orbital degrees of freedom at each site.
As the non-magnetic state has space group R 3m (No. 166), we impose time reversal T = is y K, inversion P = τ z , and threefold C 3z = exp Àiπs z =3 À Á and twofold C 2x = − is x rotational symmetries, where s i and τ i are Pauli matrices for spin and orbital, respectively. The onsite Hamiltonian satisfying all the symmetries of the nonmagnetic state is h 0 = e 0 + e 5 τ z . Along the z direction, the nearest-neighbor hopping matrices are T 4 t j + a 4 ,j = t z 0 + it z 3 s z τ x + t z 5 τ z , where a 4 = (0, 0, a z ), a z is the out-of-plane lattice parameter. For the in-plane directions, the hopping matrices are t j + a 1 ,j = t 0 + it 1 s x τ x + it 4 τ y + t 5 τ z = ðt j,j + a 1 Þ y , t j + a 2 ,j = C 3z T 1 C À1 3z , and t j + a 3 ,j = C 3z T 2 C À1 3z , where a 1 = (a, 0, 0), a 2 = C 3z a 1 , a 3 = C 3z a 2 , a is the in-plane lattice parameter (see Methods for further details). We consider the effect of the layer-alternating (i.e., A-type) antiferromagnetism by adding (−1) n−1 mσ z to h 0 , where n is the layer index. While this term breaks time reversal symmetry, inversion symmetry remains the symmetry of the lattice system. However, finite even-layer systems break inversion symmetry and have non-zero axion angle according to Eq. (10). Figure 4a shows the orbital part of the axion angle calculated with the tight-binding parameters of MnBi 2 Te 4 derived in ref. 36 (We make a momentum-dependent overall energy shift to obtain an insulating filling at half filling as we describe in Methods. The band structure and electric susceptibility are shown in Supplementary Figs. 1 and 2.). At high energies above 1 eV, the axion angle is well approximated by the optical layer Hall conductivity for any number of layers. However, the axion angle deviates significantly from the optical layer Hall conductivity as the photon energy gets lower below 1 eV. The deviation at the low energy increases with the number of layers because the surface massive Dirac fermion is then more localized at the surfaces and increases the second term in Eq. (10), making the static axion angle reach the quantized value θ = π.
T xx and spin magneto-electric coupling are much smaller than the axion angle, as shown in Fig. 4b,c, respectively. This is consistent with our expectation that these origin-independent magneto-electric couplings are strongly suppressed in systems with local inversion symmetry. Furthermore, they decrease inversely with the number of layers N l 31 , because only O(1/N l ) portion of layers near the top and bottom generates a nontrivial response. This contrasts to the case with a finite T xx for N l → ∞, where the response is coming from O(N l ) layers.
As T xx is relatively small, optical axion electrodynamics dominates the Kerr and Faraday effects in this system. Figure 5 shows the Kerr and Faraday rotation angles calculated with the magneto-electric coupling. We consider the case where the model system (medium 2) is encapsulated by medium 1 and medium 3, having frequency-independent refractive indices n 1 = 2.2 and n 3 = 2.4, respectively, corresponding to those of the hexagonal Boron nitride and diamond at photon energy around 1 eV 42,43 . The calculated Kerr rotation angle φ K is about 0.02 ∘ at photon energies larger than 1 eV, which is about one order of magnitude smaller than φ K ≲ 1 ∘ in typical ferromagnets although our antiferromagnetic system has zero net magnetic moment. The Kerr angle in real MnBi 2 Te 4 can even be much enhanced because of the contributions from higher-energy bands that we do not include here. As we consider n 1 ≠ n 3 , the Faraday rotation is nonzero because the cancellation between the top and bottom surfaces is incomplete. The Faraday effect is two orders of magnitude weaker than the Kerr effect.

Discussion
Our theory of optical axion electrodynamics fills a crucial missing piece in the macroscopic theory of magneto-optic effects in antiferromagnets developed mostly by Graham and Raab 21,[25][26][27] . They used only origin-independent T ij and S ijk to ensure physically meaningful results such as the consistency with the reciprocal relations. However, our approach shows that we can include the origin-dependent axion magneto-electric coupling in the theory, and it is precisely the axion angle that controls the Kerr effect in antiferromagnets with local inversion symmetry. As we show by using the low-energy tight-binding model of MnBi 2 Te 4 , the omission of the axion electrodynamics can underestimate the Kerr effect by orders of magnitudes. In general, the same suppression of T ij and S ijk is expected in systems with bulk symmetries that reverses an odd number of spacetime coordinates. As long as those symmetries are broken at the surfaces, the axion magneto-electric coupling is not much affected by the symmetries, such that axion electrodynamics dominates the response.
In fact, the trace part of the magneto-electric coupling was included in the study by Hornreich and Shtrikman 20 prior to refs. 21,25-27, However, their estimate of the effect was four orders of magnitudes smaller than the value experimentally observed subsequently 44 . This inconsistency lead to the appearance of theories based on different approaches [24][25][26][27] , all of which do not include the axion electrodynamics. Our introduction of the surface-sensitive magneto-electric coupling and the study of double interfaces allow for precise quantitative understanding of the Kerr effect including the axion electrodynamics, especially for thin films.  Real and imaginary part of the orbital magnetoelectric axion angle. Bulk LH indicates the optical layer Hall conductivity calculated with periodic boundary conditions. Bulk LH approaches the exact axion magnetoelectric coupling as photon energy increases. Orbital magneto-electric T xx and the spin magneto-electric coupling. Spin parts contribute to axion angle and T xx through θ ðzÞ,s = ð2G s xx + G s zz Þ=3 and T s xx = ðG s xx À G s zz Þ=3. c Spectra for two layers. These are three orders of magnitudes smaller than the orbital axion magnetoelectric coupling. d Layer-number dependence at photon energy ℏω = 2 eV. Solid and dashed lines are curves fitted with a/N l + b form, where a and b are fitting parameters. All optical response functions are calculated with γ = 10 meV to broaden the resonance through ω → ω + iγ.
While we focus on the Kerr effect due to macroscopic magnetoelectric coupling for fully compensated antiferromagnets, there is another microscopic mechanism based on the spatial modulation (phase change as well as the attenuation) of the electric field EðzÞ = E 0 e ik z z 24 , where k z = nω/c is complex valued. However, the microscopic mechanism produces only minor effects. To see this, let ϕ K,0 be the Kerr rotation angle by a single layer with net magnetization in A-type antiferromagnet. Considering the reflection of each layer as well as the complex phase rotation during the propagation, one obtains the Kerr angle ϕ K = ϕ K,0 ð1 À e 2ik z a z + e 4ik z a z À e 6ik z a z + . . .Þ= ð1 + e 2ik z a z + e 4ik z a z + e 6ik z a z + . . .Þ = À ik z a z ϕ K,0 + Oððk z a z Þ 2 Þfor an even number layers 24 , where a z is the layer spacing. This Kerr angle is smaller than the axion-induced ϕ K ≈ ϕ K,0 because a z ≪ λ for the wavelength down to the UV regime.
In the static limit, both Faraday and Kerr effects are often considered as manifestations of the axion electrodynamics 8,13,14 . It is because the systems under consideration have finite net Hall conductivity. The same sign of the Hall conductivity is induced on the top and bottom surfaces of a Z 2 topological insulator by either external magnetic fields 8,13 or coupling to ferromagnets 14 . The main focus of those studies is the manifestation of the half-quantized surface Hall conductivity, rather than the magneto-electric response of antiferromagnets.
An open question we leave for future studies is to formulate a quantum geometric theory of optical axion electrodynamics in periodic systems, generalizing the Chern-Simons integral in the static limit. A drawback of calculating intra-cell optical magneto-electric coupling through Eq. (9) (or calculating the layer Hall conductivity) is that it does not capture the topological magneto-electric effect because we drop the second term in Eq. (10). A unified formula that captures both intra-cell optical magneto-electric coupling and topological magnetoelectric effect is desired, and it is likely to require extending the quantum geometric formulation for electric dipole moments 45 to magnetic dipole and electric quadrupole moments and defining a proper optical Chern-Simons integral. However, this may not be feasible, in which case we are forced to work with quasi-two-dimensional systems.
Finally, we note that the optical axion angle we define should be distinguished from the dymanical axion fields 46 . In our optical axion electrodynamics based on linear response theory, the effective action S OA ∝ ∫ dωdxθ(ω, x)E(ω, x) ⋅ B(ω, x) for non-absorptive media describes the propagation of the electromagnetic fields modified by elastic scattering by the medium. The optical axion angle is a ground-state property, which is non-dynamical. On the other hand, the dynamical axion field interacts spacetime-locally with the electromagnetic fields through S dynamic ∝ ∫ dtdxθ(t, x)E(t, x) ⋅ B(t, x). This describes Raman scattering where the energy or momentum of the incoming electromagnetic field is tranfered to the dynamical axion field. The interplay between the two distinct phenomena is an interesting research direction.

Generalization to include natural optical activity
Electromagnetic multipole moments. Let us consider electric dipole Pi , electric quadrupole Q ij , and magnetic dipole M i moment densities induced by electric and magnetic fields 21 : for monochromatic electromagnetic fields in time domain, where P 0 i , Q 0 ij , and M 0 i are permanent multipole moments, a ijk = a kij , a 0 ijk = À a 0 kij , G ij = G ji , G 0 ij = À G 0 ji , The ellipsis "…" indicates electric-octupole/magnetic-quadrupole or higher-order multipole contributions that we neglect here. Here, the primed susceptibility tensors transform oppositely under time reversal compared to the non-primed ones. For example, while χ ij is even under time reversal, χ 0 ij is odd under time reversal.
Electromagnetic multipole moment densities are defined by 21 where we split the orbital magnetic and spin parts, which respectively originates from the minimal coupling ∇ → ∇ + ieA and the explicit dependence on B independent of the minimal coupling.
The origin dependence shows the ambiguity of defining the surface degrees of freedom. Let us recall that we define the surface current density j s = R d s =2 Àd s =2 dzj i ðzÞ over the region − d s /2 ≤ z ≤ d s /2. While we can define a physically meaningful value of d s based on the surface property of a system, this value is still not completely uniquely defined (e.g., if d s is the surface thickness, 1.01d s also makes sense as a surface thickness). Because of the ambiguity of d s , the amount of the bulk-conductivity contribution to j s can also vary. For example, let us consider the xy component of the surface conductivity. It has a bulk-conductivity term as well as the magnetic dipole and electric quadrupole contributions. Let us suppose that the surface is defined as the interface between medium 1 (at z < 0) and medium 2 at (z > 0).
where δG ðzÞ is the change ofG ðzÞ yy by the shifting of the origin by d = (0, 0, d z ). This shows that, while we can define σ s xy ðωÞ as the difference ofG ðzÞ yy between two media for any value of d s , we have to shift the spatial origin by − d s / 2 and d s /2 for medium 1 and medium 2, respectively.

Quantum mechanical expressions of the magneto-electric coupling
Linear response theory. The susceptibility tensor for the linear response of an operatorÂ to the external field F B AðtÞ =ÂðtÞ D E is given byχ whereB = ∂ĤðFÞ=∂F B | F B = 0 is conjugate to the external field. Here, we take the Heisenberg pictureÔðtÞ = e iĤtÔ e ÀiĤt . In the frequency domain, the susceptibility tensor can be written asχ where |ni is the energy eigenstate of the unperturbed Hamiltonian with energy E n = ℏω n , ω mn = ω m − ω n , and f nm = f n − f m is the difference between the Fermi-Dirac distribution function f of the |ni and |mi states, and the FS terms originate from the Fermi surface, i.e., they have momentum-space derivatives acting on the Fermi-Dirac distribution function (∂ k f n ) in the momentum space representation. The derivation goes as follows. At zero temperature, we havẽ dt 0 e iωðtÀt 0 Þ Θðt À t 0 Þh½ÂðtÞ,Bðt 0 Þi dt 0 e iðω + iΓÞðtÀt 0 Þ f n ðhn|ÂðtÞ|mihm|Bðt 0 Þ|ni À hn|Bðt 0 Þ|mihm|ÂðtÞ|niÞ dt 0 e iðω + iΓÞðtÀt 0 Þ f n ðe Àiω mn ðtÀt 0 Þ hn|Â|mihm|B|ni À e iω mn ðtÀt 0 Þ hn|B|mihm|Â|niÞ where we introduce a finite relaxation rate Γ for convergence of time integral,Ô n =P nÔPn withP n = |ni n h | is the projection ofÔ to states |ni, and P n f n hn|½B n ,Â n |ni = P n≠m f nm A nm B mn . It is often convenient to separate the real and imaginary parts of hn|Â|mihm|B|ni as follows.
where we assumeÂ andB are Hermitian operators. Let us take electric susceptibilityχ ij as an example, which is defined by Then,Â =P i is the electric polarization density andB = ÀP j V is the polarization, so we havẽ Similarly, other susceptibility tensors are given bỹ n,m f nm ω mn ωðω mn À ωÞ hn|Q ij |mihm|P k |ni = a kij + ia 0 kij , Structure of the optical magneto-electric coupling. The bare magneto-electric coupling can be decomposed into three parts as follows.
The K-term and C-term correspond to the Kubo-like and Chern-Simons terms in ref. 19 of the static limit. The last A-term is nonzero only when ω ≠ 0. While the C-term was called the Chern-Simons term, it has additional terms, actually. In the static limit, Here, θ CS is the axion angle given by the Chern-Simons integral in the three-dimensional Brillouin zone 17,18 θ CS = À 4π 2 3V ϵ ijk ImTr Pr i Pr j Pr where P is the projection to the occupied states, and A k n 0 n = hu n 0 |i∂ k |u n i is the non-abelian Berry connection for the occupied states. Since the last term vanishes in insulators with vanishing Chern number, only the Chern-Simons integral remains in the expression. The quadrupole term takes the form lim ω!0 ωa 0 kli ðωÞ Â Ã = À e 2 _V ϵ jkl ImTr½Pr k Qr lri . Lastly, g ij = P n2occ,m2unocc r i nm r j mn is the quantum metric of the occupied states. The quantum metric term in Eq. (41) cancels the Fermi surface contribution from the quadrupole term (the quantum metric contribution in electric quadrupole responses was discussed in 47,48 ), such that there is no Fermi surface contribution to the magnetoelectric coupling G ij . In comparison, note that its time-reversalsymmetric counterpart G 0 ij has a Fermi surface contribution, which gives rise to natural optical activity in metals, termed gyrotropic magnetic effect 49,50 . The quadrupole and quantum metric terms were missed in previous studies 17,18 , but they do not affect the axion angle.
Let us derive Eq. (42). A key equation in our derivation is which follows from the Wannier representation: Using this identity, we can show that It follows that This is equivalent to Eq. (42). At finite ω, the trace part of the magnetoelectric coupling is not represented as a Chern-Simons integral by the same approach.
Gauge-invariant expressions. SinceG Here we derive such gauge-invariant expressions. We focus on the orbital magneto-electric coupling to simplify expressions. It is straightforward to include spin parts as the spin magnetic moment is well defined in momentum space. The expressions of T ij , T 0 ij and S ijk are also given for completeness.
Decomposition of the position operator Here we derive Eq. (9). Let us consider the matrix element of the position operator in the Bloch state basis. We transform it to the Wannier basis |w αR by where we use |w αR i =T R |w α0 i in the second line, whereT R is the translation operator by R, and useT y R 0r zT R 0 =r z + R 0z and the orthonormality of the Wannier states hw βR 0 |w αR i = δ αβ δ R,R 0 in the third line.
To get Eq. (9) from Eq. (57), we define A z βα ðkÞ = P R hw β0 |r z |w αR ie ikÁR , such that where N is the number of k points. Then, the first term in the last line of Eq. (57) becomes and we arrive at Eq. (9).
The second term is nonzero only when n = m because where we use that ∂ k z δ k,k 0 ≠ 0 requires k 0 ! k, and P α ψ mk |ψ αk ψ αk |ψ nk = δ mn .

Macroscopic electrodynamics in the medium
Electric displacement D and magnetic field H satisfying Maxwell's equations whereF = F À iF 0 for F = A, T, U, X, and up to electric quadrupole-magnetic dipole order. While Maxwell's equations do not uniquely specify the form of D and H, additional requirements from the reciprocity relations and spatial-origin independence gives Eq. (64) as shown in 25 . Here, the free charge and current ρ f and J f are boundary charge and current appearing due to the change of material properties across the interface that is not described byT ij ,Ũ ij , and S ijk = ða 0 ijk + a 0 jki + a 0 kij Þ=3. In PT-symmetric systems where G 0 ij = a ijk = 0, Wave equation. The wave equation up to electric quadrupole/magnetic dipole takes the following form 21 : whereχ ij = χ ij À iχ 0 ij satisfying χ ij = χ ji and χ 0 ij = À χ 0 ji , κ i = k i /|k| is the propagation direction of light, n is the refractive index,σ ijk = σ ijk À iσ 0 ijk is the complex bulk conductivity coefficient defined bỹ σ ij ðqÞ =σ ij ð0Þ +σ ijk q k + . . ., and Let us assume C 3z and C 2x symmetries for simplicity. We further impose that the bulk Hall response is zero, i.e., χ 0 ij = 0, in order to focus on the magneto-electric and electric-quadrupole effects. For κ = ±ẑ, the wave equation is The refractive index satisfying the wave equation is given by for circular polarization± =x + iŷ.
Reflection and transmission from a single interface. We consider the interface of medium 1 (z > 0) and medium 2 (z < 0) with the surface normalẑ. For normal incidence, in the circularly polarized basis, within the media μ = 1 or 2, where n ± depends on the sign κ z , and κ z = − 1 for incident and transmitted light, while κ z = 1 for reflected light. Here, As we consider light incident from medium 1 to medium 2, the electric field in medium 1 consists of incident and reflected fields while that in medium 2 is the transmitted field.
where r = r + + r + + by C 3z symmetry and the continuity of E at the interface.

Data availability
The data that support the findings of this study are available in the main text and Supplementary Information. Further information is available from the corresponding author upon reasonable request.